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partition function. We also derive a formula for the exchange probability in the replica exchange 
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I. INTRODUCTION 

One of the most challenging subjects in the field of 
computational physics is to develop efficient methods for 
long-range interacting systems. The difficulty of long- 
range interacting systems comes from the fact that the 
number of interactions rapidly increases with increasing 
the number of elements TV of the system. For exam- 
ple, in the case of systems with pairwise interactions, the 
number of interactions is proportional to A^ 2 . Therefore, 
if one uses a naive simulation method in such systems, 
the computational time per one step rapidly increases 
in proportion to A^ 2 , which is quite contrast to the case 
of short-range interacting systems in which the compu- 
tational time increases in proportion to N . In order to 
overcome the difficulty, many simulation methods have 
been proposed until now [lHlOj. 

Recently, the author and Matsubara have developed 
an efficient Monte Carlo (MC) method called Stochas- 
tic CutOff (SCO) method for long-range interacting sys- 
tems [Til] - The basic idea of the method is to switch 
long-range interactions Vij stochastically to either zero 
or a pseudointeraction Vij with the detailed balance con- 
dition satisfied. Interactions are switched by u sing the 
Stochastic Potential Switching (SPS) algorithm [ll 
Because most of the distant and weak interactions are 
eliminated by being switched to zero, the SCO method 
greatly reduces the number of interactions and computa- 
tional time in long-range interacting systems 
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thermore, this method does not involve any approxima- 
tion because the detailed balance condition is satisfied 
strictly. Fukui and Todo have recently developed an ef- 
ficient MC method based on similar strategy by use of 
different pscudointeractions and different way of switch- 
ing interactions fl5j ] . 

In the present work, we reformulate the SPS algorithm 
which is used in the SCO method. This reformulation 
gives us a generalized version of the Fourtuin-Kasteleyn 
representation of the partition function in the Ising fer- 
romagnetic model [l6l . [l7| . This representation of the 
partition function is used to derive new formulae for in- 



ternal energy and heat capacity measurements. Since 
these formulae consist of only terms which survive as Vij 
in the above-mentioned potential switching process, the 
computational time for the measurements are reduced 
greatly. We also derive an formula which reduces the 
computational time to estimate exchange probability in 
the replica exchange MC method [l8|. In order to ver- 
ify the formulae, we perform some MC simulations in 
three dimensional magnetic dipolar systems. The results 
clearly show the validity of the formulae. 

The organization of the paper is as follows. In £ fl"fl we 
reformulate the stochastic potential switching algorithm. 
A generalized Fourtuin-Kasteleyn representation of the 
partition function is presented in this section. In ijllll 
and i]IVl we derive formulae for internal energy and heat 
capacity measurements and a formula for the replica ex- 
change MC method, respectively. The validity of these 
formulae is confirmed numerically in WII Section [VIII is 
devoted to a summary and discussions. 



II. REFORMULATION OF THE SPS 
ALGORITHM 

Before reformulating the SPS algorithm, we briefly ex- 
plain the SPS algorithm [l^, EH . We hereafter consider 
a system with pairwise long-range interactions described 
by the Hamiltonian 



H 



E 



KJ 



Vij (Si, Sj), 



(1) 



where Si is a variable associated with the i-th element 
of the system. In this algorithm, Vij is stochastically 
switched to either Vij or Vij with a probability of or 
1 — Pij , respectively. The probability P,j is 

PijiSi, Sj) = exp[/3(A^-(S.„ Sj) - AV*)], (2) 

where /3 is the inverse temperature, 

AVij IS, ,Sj) = Vij{Si ,S 3 )- Vij{Si ,Sj), (3) 
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and AV^* is a constant equal to (or greater than) the 
maximum value of AV^j(Sj, Sj) over all Si and Sj. We 
can choose the potential Vij arbitrarily. On the other 
hand, using Py(S^, Sj), the potential Vij is given as 

VaiSi, Sj) = VijiSi, Sj) - /T 1 log[l - Py(5i, 5,-)]. (4) 

With this potential switching process, the algorithm pro- 
ceeds as follows: 

(A) Potentials Vij are switched to either Vij or Vij with 
the probability of Pjj or 1 — Pjj , respectively. 

(B) A standard MC simulation is performed with the 
switched Hamiltonian %' expressed as 



n> = E' ^ s i) + E" ^ 



(5) 



where runs over all the potentials switched to 
V - and X)" runs over those switched to V. The 
potential is fixed during the simulation. 

(C) Return to (A). 

It is shown that this MC procedure strictly satisfies the 
detailed balance condition with respect to the original 
Hamiltonian of eq. (JTJ. 

In the SCO method, Vij is set to zero to reduce the 
computational time of %' defined by eq. (0). Further- 
more, the use of an efficient method enables us to reduce 
the computational time of the potential switching proce- 
dure (A) greatly (see ref. [IH for details). 

In order to reformulate the SPS algorithm, we first 
introduce a graph variable gij defined by 



if Vij is switched to Vij , 

1 if Vij is switched to Vij . 



(6) 



We next introduce a weight u)(Si, Sj - ,gij) defined by 



w ij {Si , Sj ; g%j) 



(7) 

This weight is analogous to the weight introduced by Ed- 
wards and Sokal [19j]- We hereafter show that the SPS 
algorithm is a MC method which realizes equilibrium dis- 
tribution defined by 



Psps{{S I },{g i] }) = Z S p S Y[ .^.Uij{Si,Sj-,9ij), (8) 
where 



Zsps = Tr 



II w «( 5i ' S i>9ij)- (9) 



As it is readily derived from eqs. (J2J), (J3J), (jlj), and (0, 
the weight u>ij{Si, Sj] gij) has the following property: 



Tr^—ciWy^Si, Sj\gij) 



eM-/3Vij{Si,Sj)}. (10) 



This equation naturally leads us to the following new 
representation of the partition function 

Z{p) = Z S ps(^) = Yv :S . . ;/ } \[^ ^,,iS,.S,:g,,). 

(11) 

We also find that the probability that some configuration 
{Si} is realized in the SPS algorithm is given as 

P{{Si}) = Tr {ffi . } PsPs({Sy,{<fe}) = P B {{Si}), (12) 
where Pb is the Boltzmann distribution defined as 

P B ({5J) = Z(/3)- 1 exp[-/3^. y ij {Si,S j )\. (13) 

L 1 — <i<3 J 

This is the reason why the Boltzmann sampling is 
achieved by the SPS algorithm. 

We next show that the equilibrium distribution of the 
SPS algorithm is given by eq. flSJ). Let us denote the 
transition probability in the step (A) of the SPS algo- 
rithm as W\{{gij} — > {ffijlK'S'i}) anci lh at m tne s t e P 
(B) as W B {{S l } -> {S'i}\{gij}). It should be noted that 
the procedure in the step (A) updates the graph vari- 
ables {gij} with fixing the configuration variables {Si}, 
and that in the step (B) updates {Si} with fixing {gij}. 
In the following, we will show that the two transition 
probabilities satisfy the detailed balance conditions 

Psps({^},{^})Wa({^} Wij}\{Si}) 

= P S Ps({Sy, Wij})W A {{g'ij} -> {QiMSi}), (14) 

P SPS {{S t },{g 13 })W B {{S t } H- {S'^gij}) 

= Psps({^}, {g l3 })W B {{S' l } -> {S^gij}). (15) 

These two equations clearly show that the equilibrium 
distribution of the SPS algorithm is Psps{{Si} , {g^}) . 

We start from the proof of eq. (| 14|) . It can be easily 
seen from the procedure in the step (A) that 

W A {{9ij} -+ Wij}\{Si}) 

= n (0) ' Paw, Sj) n (1) ' [i - Pi^Sj)], (i6) 

where the product Y[ runs over all the pairs with g[j = 

and n' 1 '' runs over those with g[j = 1. It should be 
noted that W A does not depend on {gij}- To rewrite 
the right hand side of the above equation, we utilize the 
following two equations: Firstly, it is found from eqs. ([2"]l. 
©, and that 

Pij{Si, Sj) = u>ij{Si, Sj^g'ij = 0)eP v ^ S *' S i\ (17) 
Secondly, we can rewrite 1 — Py (Sj, Sj) as 

1 — Pij {Si , Sj ) 

= {e-^ S - S ^ -UijiSuSj^j = Q))^S i ,S i ) 
= u k j{S i ,Sj;g' i j = iy v ^ S *> S J\ (18) 
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where we have used eqs. (| 10[) and (fl7|) . By substituting 
these two equations into eq. (|16p . we obtain 

W A {{9ij} -> MM 8 *}) 

= ZsrsPsrsttS^^g'^U e? v ^ 8 ^. (19) 

This equation shows that W-A satisfies the detailed bal- 
ance condition (|14[) . 

Our second task is to prove eq. (p~5]) . Since we perform 
a MC simulation with the switched Hamiltonian H', the 
transition probability Wb satisfies the detailed balance 
condition 

e -/8{E w ««(5 4 ,5 i )+E <1) * S„S | 

x^ B ({^}^{5'J|{. 9y }) 
= e -^{E <0) ^-(^.S^+E' 1 ' v. A" A" } 

xWb({^} -+ {Si}|{<fe}), (20) 

where the sum J^ ^ runs over all the pairs with = 
and y~] runs over those with gij = 1. By multiplying 
e -/3E ( ' Av ij to the both hand sides of the equation, we 
obtain 

WB({Si} {s'Mgi^Ui^^^m) 

where we have used eq. ([7]). It is clear from eq. (J5) that 
this equation is equivalent to the detailed balance condi- 
tion (fT5|) . 

We next turn to the new representation of the parti- 
tion function, i.e., eq. (jlip . This is a generalization of the 
Fourtuin-Kasteleyn representation of the partition func- 
tion in the Ising ferromagnetic model [HI, Hill- ^ n f ac t> ^ 
is shown in the appendix A that the original Fourtuin- 
Kasteleyn representation is derived from cq. ((lip . This 
representation is more comprehensive than the original 
one in the following two points: 

1) In the new representation, potential Vij can be cho- 
sen arbitrarily. On the other hand, the original 
Fourtuin-Kasteleyn representation corresponds to 
a special case in which Vj is zero. 

2) The new representation is valid no matter whether 
the configuration variables {Si} are discrete or con- 
tinuous. This is contrast to the original represen- 
tation which is derived for systems with discrete 
variables. 

This generalized Fourtuin-Kasteleyn representation will 
be used in the next section to derive formulae for internal 
energy and heat capacity measurements. 



III. FORMULAE FOR INTERNAL ENERGY 
AND HEAT CAPACITY 

In order to derive formulae for internal energy (E) and 
heat capacity (C), we use the following thermodynamic 
relations: 



(E) = -Z 



(C) = k B P 2 Z~ 



d 2 z 



_ 1 dZ_ 

9/3' 



(22) 



(23) 



As shown in the Appendix B, the formulae are obtained 
by substituting our generalized Fourtuin-Kasteleyn rep- 
resentation eq. (fTTj) into these equations and calculating 
the derivatives. The results are 



(E) = E c 



/ SPS 



(24) 



where 



(25) 



(O)sps = Tr {gi . }){Si} 0Psps({Si},{<fe}), (26) 



E const -/Z k<l AV kl' 



(27) 



i-p, 



d/3 ~ ^ 



l.'<l 



AV kl - AV, 



hi 



1 - Pm 



(29) 



Pki in eqs. (|28j) and (|29|) is the switching probability de- 
fined by cq. ((2|). It is quite important to note that the 
average in MC simulations is equivalent to (0)sps 5 



(O)mc = (O) 



SPS- 



(30) 



This comes from the fact that MC simulation with 
the SPS algorithm samples states with the probability 

In the case of Vj = 0, E const , E, and M in eqs. ([27]). 
and (f2l?|) are reduced to the following forms: 



V k i - V, 



E=y 5 ah , i , 

^k<i a "'^ V \-p u 



kl 



(31) 



(32) 
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k<l 



1 - Pm 



p, 



(33) 



where V*j is a constant equal to (or greater than) the 
maximum value of Vij (Si, Sj) over all Si and Sj. The 
point of the formulae is the presence of 5 gklt i in E and 
In general, the strength of pairwise long-range in- 
teractions Vij decreases with increasing the distance . 
Therefore, g^ for most of distant and weak interactions 
becomes zero (see eqs. ^ and &§)■ This means that 
the computational time of E and its derivative is much 
shorter than that of X)fc<; Vki which is needed in naive 
internal energy and heat capacity measurements. Al- 
though the calculation of the constant E const requires 
0(N 2 ) computational time, it is enough to calculate it 
just once at the beginning of the simulation. 

In summary, the use of the formulae (eqs. (|24l) . ([2~5"T) . 
(|3"Tj) . (f3"2")) . and (|33[0 enables us to reduce the computa- 
tional time for internal energy and heat capacity mea- 
suremen ts g reatly. This method can be used in the SCO 
method [ll| because Vij is set to be zero in the method. 



IV. A FORMULA FOR THE REPLICA 
EXCHANGE MC METHOD 

We first explain the replica exchange MC method [l8| 
briefly. This method is quite efficient for systems with 
rugged energy landscape such as spin glasses. In this 
method, we consider a system with M independent repli- 
cas and M different temperatures. The M replicas have a 
common Hamiltonian H(S). The purpose of the method 
is to sample states of the M replica system with the fol- 
lowing equilibrium distribution 

nM 
P B (S m ,p m ), (34) 
m— 1 

where S m denotes the set of real variables {Si} of the 
m-th replica and Pb is the Boltzmann distribution de- 
fined by eq. (Ti~3"|) . In the replica exchange MC method, 
the equilibration is accelerated by exchanging the replica 
at temperature T m for that at T n . By employing the 
Metropolis method, the probability of accepting the ex- 
change is given as 



W B (S, p m ; S', p n -> S', p m ; S, p n ) 

Qb(- • • ; S', /3 m ; • • • ; S,p n ; • • • ) 



= min < 1 



where 



Qb(- • • ; S,/3 m ; • • • ; S ,/?„; • • • ) 
= min[l, A B ], 



X B = exp[(/3 m - p n )(H(S) - H(S'))]. 



(35) 



(36) 



A problem of the replica exchange MC method in long- 
range interacting systems is that it costs 0(N 2 ) compu- 
tational time to calculate the exchange probability since 



Xb in eq. (|36|) contains T-L which consists of 0(N 2 ) pair- 
wise interactions. In order to overcome the difficulty, we 
consider a replica exchange MC method whose equilib- 
rium distribution is 



Qsps(£i,<7i,/3i; • • • ; Sm,9m, Pm) 

nM 
. PsPs(Sm,9 m , Pm), 
m.— \ 



(37) 



where g m denotes the set of graph variables {gij } of the 
m-th replica and Psps is defined by eq. (JSj). It should 
be noted that this replica exchange MC method sam- 
ples configuration (Si, • ■ ■ , Sm) according to the original 
probability Qb since Qsps is related to Qb as 



Tvg it ... t g M QSFS(SU91,PV,- 

= Qb(Si,Pi; • ■ ■ ; Sm,Pm), 



(38) 



where we have used eq. (fT2"| . The probability of accepting 
the replica exchange is given as 

W SPS (S,g,p m] S',g',p n -> S',g',p m ;S,g,p n ) 

. Qsps(- ■ ■ ; S', g', p m ; ■ ■ ■ ; S,g, p n ; ■ ■ ■ ) 
mm < 1 



' Csps("- ;S,g,p m ;--- ;S',g',p n ;---), 
min[l,X SPS ], (39) 



where 
Asps : 



TT (0) ' 

exp 



{Pn-PmHVijiS'^S'A + AV*,} 



x [] (0) exp [{p m - P n }{%{Si, Sj) + AV*} 
x I]' 1 " exp [{Pn ~ P m }V tl (S>, 5;.)] t"p J( 5^'^ 

1 — ~ij K°ii °jl Pn) 

x [{Pm ~ PnW^Si, Sj)} IZ^^sM ' 

(40) 

In this equation, the product runs over all the pairs 

with gij = a and J^"' runs over those with g'^ = a. 

When Vij = 0, we can rewrite the first product in the 
right hand side of eq. (|4"0)) as 



r(0)' 



J exp [{p n - p m }V*j] 



(i)' 



J exp [{p n - p m }V£] I exp [-{p n - p m }V*] 



(41) 



The second product can be rewritten in a similar way. 
By substituting these results into eq. (|4T)]) . we find 



Asps 



1 Pij(S^, S j, P m ) 
P^S'^Pn) 



JJ UJ e [{fln-/s»}{v„(5;,5J)-vg}]^ 
JJ (1) e [{0 m -M{v ij (S i ,S j )-v* i }] 1 - Pjj(Si, Sj,p n ) 



(42) 
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Since Xsps in this formula is calculated only from the 
pairs whose graph variables are one, it can be calculated 
with very short computational time as E' and in 
eqs. ([32]) and ([33]) are. 



V. THE CASE WHEN LONG-RANGE 
INTERACTIONS AND SHORT-RANGE 
INTERACTIONS COEXIST 

When long-range interactions {V^ } and short-range 

interactions {V^p} coexist, we do not need to use the 
SCO method for short-range interactions because it does 
not reduce the computational time significantly. In other 
words, we can switch Vit to or Vij with Vjj un- 
changed. This can be realized by setting = and 



(S) 



(S) 

. In this case, we can decompose E c 



§f , and X SPS in eqs. (2ft, @SJ, ((29j) , and (@DJ) as 



-^const — -^const "r" ^const> 
E = E {L} + £ (S) , 

dE _ dE^ dE^ 
~dp ~ 



(S) 



dp 



dp 



X S ps = X&l x X (s) 



L SPS 



SPS' 



E, 

(43a) 
(43b) 

(43c) 

(43d) 



where the first terms and the second terms in the right 
hand sides denote contributions from long-range inter- 
actions and those from short-range interactions, respec- 
tively. As we have already mentioned, the first terms are 
given by eqs. ([21]), ©, ([55)1. and (gSJ, respectively. On 



the other hand, by substituting Vki 
into eqs. ([27)1, ([55)1. ([2? 



V, 



(S) 



A7 



and AV^ = 



£ (s) = 

const ^? 



^ (S) = -E (fei) ^ 

9/3 



(S) 



X SPS = n <fci) CX P[^'" - PnWuiSk, S t 



X 



(S) 



and ([40)) . we readily obtain 

(44a) 
(44b) 

(44c) 
(44d) 



where we have used the fact that all the potentials are 
switched to Vm ■ Note that Pki = 1 when Vki — Vki and 
AT/* = 0. These results are quite natural because they 
coincide with the results in the usual MC procedure. 



VI. NUMERICAL TESTS 

A. Internal energy and heat capacity 
measurements 

In order to check the validity of the formulae for inter- 
nal energy and heat capacity measurements, we perform 



LU 



O 




0.2 



0.4 



0.6 



0.8 
T/D 



1.2 



1.4 



FIG. 1: (Color online) Temperature dependence of (I) inter- 
nal energy E/N and (II) heat capacity C/N in the purely 
dipolar system. Simulated annealing method is used for the 
measurement. The size N is 10 3 . Measurements are done in 
the three cases (see text). 



MC simulations of a three dimensional magnetic dipolar 
system on a L 3 simple cubic lattice. The boundary condi- 
tion is open. The Hamiltonian of the system is described 
as 



H = H 



dip 



Sj ' Sj — 3(Sj ■ Tij)(Sj ■ Vi 



(45) 

where Si is a classical Heiscnbcrg spin of |Sj| = 1, r^- is 
the vector spanned from a site i to j in the unit of the 
lattice constant a, nj = |ry|, and D = (g/j,&S)/a 3 . We 
hereafter call the system purely dipolar system. 

In simulations with the SCO method, we regard each 
term in the right hand side of eq. (|4"5]) as Vy. Since 
% = in the SCO method, AV# defined by eq. © 
is equal to 

2D/rfj when S, and Sj are anti-parallel along r i3 -. We 
therefore set AV^, which should be equal to or greater 
than AVa over all Sj and Sj, to 2D/r 3 j. By substituting 



Interaction Vij has the maximum value 



these into eq. ([2j, we obtain 



Pij = exp 



f3D | Si-Sj-m ■ rij )(Sj ■ riJ ) - 2 ' 



' ij 



(46) 
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T/D 

FIG. 2: (Color online) Temperature dependence of the av- 
erage number n of potentials per site which are switched to 
%• 

Pseudointeraction Vij is given by substituting Pij of the 
above equation into eq. (U]). 

To make a comparison between the new measurement 
methods and the usual ones, we do simulation in the 
following three cases: 

(a) Usual MC method with usual measurements. 

(b) SCO method with usual measurements. 

(c) SCO method with measurements by using the for- 
mulae derived in §IIII 

In all the three cases, simulated annealing method is 
used. The system is gradually cooled from an initial tem- 
perature T = l.ZD to 0.3L> in steps of AT = 0.05D. The 
system is kept at each temperature for 2 x 10 6 MC steps. 
The first 10 6 MC steps are for equilibration and the fol- 
lowing 10 6 MC steps are for measurement. The size N is 
10 3 . The size is not so large because we do simulation not 
only with the SCO method but also with a usual method 
which requires 0(N 2 ) computational time per one MC 
step. 

In Fig. [TJ internal energy E and heat capacity C mea- 
sured in the three cases are plotted as a function of tem- 
perature. We see that all the data nicely collapse into 
a single curve, showing the validity of the formulae de- 
rived in i jllll We also see that heat capacity has a peak 
around T/D = 0.55. This result is consistent with pre- 
vious works which show the existence of a phase transi- 
tion around this temperature [2(], [HJ . Figure [2] shows 
the average number n of potentials per site that sur- 
vive as Vij in the potential switching process. Though 
n increases with decreasing temperature, n k 70 even 
at the lowest temperature. This means that more than 
90 percent of interactions are cut off by being switched 
to Vij = 0. It is worth pointing out that the SCO 
method becomes more efficient with increasing the size. 
In fact, in the study of two-dimensional magnetic dipo- 
lar system with dipolar interactions and ferromagnetic 




50 100 150 200 



t 

FIG. 3: (Color online) Time autocorrelation functions of the 
four observables E\j, Ctj, -En, and Cn (see eq. (|51[) for their 
definitions). The size A^ is 10 3 and the temperature T is 0.3D. 
All the four observables are measured with the SCO method. 



exchange interactions , it has been found that the in- 
crease of n with size is very slow and n is 22.5 even when 
N = 256 2 = 65,536. 



B. Statistical error of the new measurement 
methods 

We notice from Fig. [T] (II) that the data in the case 
(c), i.e., those which are measured with eq. (|25l) . fluc- 
tuate more than the other data. To get some insights 
of this behavior, we consider estimating statistical er- 
ror of observables. We suppose that an observable Q 
is successively measured N times in a MC simulation 
to estimate the average jj^2^ = iQn- We assume that 
the measurement is done every 5t steps. Although this 
average value is close to the thermal average value (Q), 
they are slightly different because the number of the mea- 
surement N is finite. We hereafter call the difference 
SQ = jj Ylu=i *3m — (Q) statistical error. When the pe- 
riod of the measurement NSt is much larger than the 
correlation time tq of the observable Q, the expectation 
value of the square of the statistical error ({8Q) 2 ) is ap- 
proximately evaluated as [13, H3] 

<W) 2 ) = ^[(Q 2 )-(Q) 2 ](i + 2^). (47) 

The factor (1+2^-) in the right hand side of the equation 
comes from the fact that Q^s measured successively are 
correlated with each other. From this equation, we can 
estimate the relative statistical error as 



(Q) Vn 1 
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TABLE I: Correlation time tq, relative variance, and Cse of 
the four observables Eu, Ctj, -En, and Cn (see eq. ()49p for 
the definition of Cse and eq. (|51[) for the definitions of the 
four observables). The size N is 10 3 and the temperature T 
is 0.3-D. The time interval St in eq. (|49|l is assumed to be one 
in the estimation of Cse- 







(Q 2 ) - (Q) 2 
<Q> 2 


Cse 


Q = E V 


96 


2.0 x 10~ 5 


0.062 


Q = -B N 


3.6 


5.7 x 1CT 4 


0.068 


Q = Cu 


25 


2.0 


10 


Q — Cn 





1.6 x 10 3 


40 



where 



Cse 



(Q 2 ) - (Q) 2 
(Q) 2 



1 



'St 



(49) 



In order to estimate tq in eq. (|49[) . we measure the 
normalized time autocorrelation function defined by 



0QW = 



Q(p)Q(t) 



Q 2 -{Q)' 



(50) 



where 7TT denotes the average over a sequence of the data 
obtained by a MC simulation. Since we are interested 
in how correlation times in internal energy and heat ca- 
pacity measurements are affected by their methods, we 
calculate correlation functions of the following four ob- 
servables 



(H) 



Eu(t) = H(t), 
C v {t) = k B p 2 (H{t) 
Em(t) = E(t) + E const , 

C N (t) = k B f3 2 l (E(t)-{E) 



dim ) 

dp j 



(51a) 
(51b) 
(51c) 

(51d) 



where the subscripts U and N denote the usual measure- 
ment and the new measurement by using the formulae 
derived in £1III| respectively. The results are shown in 
Fig. |3] All the four observables are measured with the 
SCO method. We see that correlation functions in new 
measurements (full symbols) are much smaller than those 
in usual measurements (open symbols). This probably 

comes from the fact that E and ^| can change without 
change in spin configurations since they depend on not 
only {Si} but also {gij}- This result shows that large 
fluctuations in heat capacity measured with the new for- 
mula do not originate from an increase in the correlation 
time. 

We next estimate the correlation time tq from </>q(£) 

as 



dt> Q (i'). 



In practice, we set the lower limit to be one and adjust the 
upper limit so that it is much larger than the correlation 
time. For example, the upper limit was set to 3, 000 when 
we estimated the correlation time for Etj . The estimated 
values are shown in Table HI We estimate the correlation 
time for Cn to be zero because it is negligibly small. We 
calculated r for Cn with changing the upper limit from 2 
to 10,000, and found that it does not exceed 0.2 for any 
upper limits. 

In Table HI we also show the values of relative variance 
and those of Cse- Wc find that the relative variances in 
new measurements are larger than those in usual mea- 
surements, especially in the heat capacity. Concerning 
the internal energy, we can explicitly show the increase 
in variance from eq. (|25p as 



k^T 2 C = (H 2 ) - (H) 1 









-<*: 




- (Es 




(52) 



(53) 



where we have used the fact that denned by eq. (f2"9")l 
is always positive. Since (En) = (Etj), the inequality 
([53]) shows that the relative variance of En is larger than 
that of En. We next turn to how Cse is affected by the 
measurement methods. We see from Tablc|I]that there is 
no significant difference between Cse for E\j and that for 
i?N- On the other hand, Cse for Cn is four times larger 
than that for Cu . This is the reason why the heat capac- 
ity in the case (c) fluctuates more than that in the other 
cases. It should be noted that the relative statistical er- 
ror is proportional to Cse (see eq. (|48|) ). Equation (|48|) 
also tells us that the number of measurements with Cn 
should be 16 times as large as that with Cu to make both 
the statistical errors the same. 

In summary, to attain a certain accuracy, estimation 
of the heat capacity with Cn requires larger number of 
measurements than that with Ctj since fluctuations in 
Cn are larger than those in Ctj. On the other hand, 
the efficiency in the internal energy measurement with 
En is almost the same as that with Etj in the present 
case. When the heat capacity measurement with Cn does 
not work well, it might be possible to estimate the heat 
capacity with lower statistical error by doing numerical 
differential of the internal energy which is evaluated with 
En- 



C. Replica exchange method 

To examine the validity of the formula derived for the 
replica exchange MC method, we again do simulations in 
the following three cases: 

(a) Usual MC method with usual replica exchange MC 
method. 
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FIG. 4: (Color online) Temperature dependence of (I) internal 
energy E/N and (II) heat capacity C/N in the purely dipolar 
system. The replica exchange MC method is used for the 
measurement. The size N is 6 3 . Measurements are done in 
the three cases (see text). 



(b) SCO method with usual replica exchange MC 
method. 

(c) SCO method with replica exchange MC method by 
using the formula derived in £|IVI 

In all the three cases, the number of MC steps for equili- 
bration and that for measurements are 10 6 . The number 
of temperatures is 40 and the size N is 6 3 . Figure 0] 
shows the result of internal energy and heat capacity 
measurements. The measurements are done with usual 
method in the cases (a) and (b), and with our formu- 
lae (eqs. (pM|) and (pp)) ) in the case (c). We again see 
that all the data nicely collapse into a single curve. This 
result clearly shows the validity of the formula derived 
in ^IVI We also see that the heat capacity in the case 
(c) fluctuates more than that in the other cases, as we 
have seen in Fig. [TJ In Fig. [5l we show the temperature 
dependence of the probability P ex (Ti) that exchange tri- 
als between i-th and i + 1-th replicas are accepted. We 
find that P cx (Ti) with our formula (case (c)) is smaller 
than that with usual method (the other cases). From this 
result, one may consider that our replica exchange MC 
method is less efficient than the usual method. However, 
this is not true because the computational time of Xgps 
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FIG. 5: (Color online) Temperature dependence of the replica 
exchange probability P ex in the purely dipolar system. The 
size N is 6 3 and the number of temperatures is 40. 
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FIG. 6: (Color online) Temperature dependence of the replica 
exchange probability P cx in the ferromagnetic dipolar system. 
The size N is 6 3 and the number of temperatures is 40. 



(cq. (|42p) is much shorter than that of X# (eq. {35}). 

We next show results of the replica exchange MC 
method when long-range interactions and short-range in- 
teractions coexist. The Hamiltonian is consist of the 
long-range magnetic dipolar interactions (eq. (|45|) ) and 
the short-range exchange interactions 



Tie 



-jy Si 



Sj (J>0), 



(54) 



where the sum runs over the nearest neighbouring pairs. 
The ratio D/J is 0.1. We hereafter call the system ferro- 
magnetic dipolar system. The SCO method is used only 
for the dipolar interactions. We again examine the three 
cases mentioned in the previous paragraph. The number 
of MC steps for equilibration and that for measurements 
are 10 6 . The number of temperatures is 40 and the size 
N is 6 3 . We first have confirmed that internal energy 
and heat capacity in the three cases coincide with each 
other. We next examine the temperature dependence of 
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TABLE II: Ergod time in the purely dipolar system and the 
ferromagnetic dipolar system. Measurements are done in the 
three cases (see text). 







purely dipolar 


ferromagnetic dipolar 


case 


(a) 


1.9 x 10 4 


3.3 x 10 3 


case 


(b) 


4.4 x 10 4 


3.5 x 10 3 


case 


(c) 


6.2 x 10 4 


4.5 x 10 3 



the replica exchange probability P cx (T; ) . Figure |6] shows 
the result. We see that the reduction of the exchange 
probability in case (c) is not as large as that of the purely 
dipolar system (Fig. [5]). This result is reasonable because 
the SCO method is used only for the long range interac- 
tions and its contribution (ATg PS in eq. (J33Ji)) is not large. 
Recall that the ratio D/J is 0.1. We also measure the 
ergod time defined by the average MC step for a specific 
replica to move from the lowest to the highest tempera- 
ture and return to the lowest one. The result is shown 
in Table [TTJ In both the systems, the ergod time in the 
cases (b) and (c) is larger than than that in the case (a) , 
meaning that the use of the SCO method increases the 
ergod time. However, the increase of the ergod time in 
the ferromagnetic dipolar system is not as large as that 
in the purely dipolar system. These results show that 
the SCO method is particularly efficient for systems with 
strong short-range interactions and weak long-range in- 
teractions. This feature of the SCO method has already 
been pointed out in the previous work [Tl| . 



VII. SUMMARY 

In the present work, we have derived useful formulae 
for the SCO method to estimate internal energy, heat 
capacity, and replica exchange probability in the replica 
exchange MC method. We can reduce the computational 
time of these quantities greatly by using the formulae be- 
cause they only contain terms which are not cut off by the 
SCO method. On the other hand, we have found that the 
use of the formulae could cause a decline in the efficiency 
of the measurement and that in the exchange probability. 
When the new methods do not work well, the analyses 
done in the present paper, such as the estimations of the 
statistical error, the replica exchange probability, and the 
ergod time, might be helpful to figure out the reason and 
to get rid of it. Anyway, we hope that these formulae 
make the SCO method more useful and attractive. 

The other achievement of the present work is the 
derivation of the new Fourtuin-Kasteleyn representation 
of the partition function, i.e., eq. (|11[) . This representa- 
tion is more comprehensive than the original one because 
our representation includes arbitrariness in the choice of 
pscudointcractions {Vy}. Furthermore, this representa- 
tion can be used no matter whether the variables {Si} 
are discrete or continuous. We hope that this representa- 



tion becomes the basis of new algorithms as the original 
Fourtuin-Kasteleyn representation lead to the Swendsen- 
Wang cluster algorithm (2JI in the Ising ferromagnetic 
model. 
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Appendix A: Derivation of the Fourtuin-Kasteleyn 
representation of the partition function in the Ising 
ferromagnetic model from eq. 

In this appendix, we consider the Ising ferromagnetic 
model whose Hamiltonian is described as 



ij) 



(Al) 



where <n = ±1 and J > 0. We set Vij = — JoiOj and 
consider a special case that Vij = and V*j = J. Then, 
we find 



Pi = \ Cxp ( -2 ^ J ) = °J')> (A2) 



and 



J - r 1 Iog[l - cxp(-2/?J)], (a % = a j ) , 

OO (<7j = -CTj). 



By substituting this equation into eq. ([7]), we obtain 

uJij((Ti,aj;gij) 



(A3) 



exp(-/3J) 
exp(fij)[l - exp(-2j8J)]5 (7 . 



(9ij 



0) , 

1) . 



(A4) 



It is important to note that the product Y[(ij) u ij * s P r0 ~ 
portional to the joint distribution of {Si} and {gij} in- 
troduced by Edwards and Sokal [l9j. In this sense, we 
can regard the product of w, 3 defined by eq. ([7]) as a 
generalization of Edwards and Sokal's joint distribution. 

As it has been pointed out in rcf. [19) . we can eas- 
ily obtain the Fourtuin-Kasteleyn representation of the 
partition function from this joint distribution. We first 
rewrite ]l<y> Wi i as 



n 



ij) 



{exp(-2/3J)} 



:{l-exp(-2/3J)}^n! 1 l^^' ( A5 ) 



ij) 



where N int = J2{tj) and = <Wi- Thc product 

Yiuj\ runs over all the pairs with <7y = 1. By using 



10 



cq. (fTTj) . wc obtain 



Z(/3) 



^N int /3J 



Tr {gi . } {exp(-2/?J)} 



:{1 - cxp(-2/?J)} A ' b Tr {ffi} fl^ <W (A6) 



Wc next consider calculating Trj CT .} ^i,^ in the 
above equation. We hereafter call pairs with <7,j = 1 
bonds and a set of sites which are connected by bonds a 
cluster. Because of the presence of <5 CTijCTj , the values of 
Cj in a cluster are forced to be the same. Therefore, we 
find 



Tr 



{Si} 



n 



(i) 



lJV c ; 



rfffy} 



(A7) 



where AT c i ustor is the number of clusters. By substitut- 
ing this equation into eq. (|A6[) . we obtain the Fourtuin- 
Kasteleyn representation of the partition function in the 
Ising ferromagnetic model, i.e., 



Tr {ffij} {exp(-2/3J)} 



:{l-cxp(-2/?J)} iVb x 2 



N int -N b 



r{9ij} 



(A8) 



where we have used the identity 



(B5) 

to go from the third line of eq. (|B4[) to the fourth. By 
substituting 



d =p-Hog[l-P kl } + f3-i T ^L 



9/3 



(B6) 



and S 9klt Q = 1 — S gklt i into eq. (|B4[) . we obtain 
From this equation and eq. (|B2[) . we find 



(B7) 



Appendix B: Derivation of eqs. (f24|) - ([29]) 



E' — —E const — E, 



(B8) 



We first derive a formula for internal energy (E). We 
see from eqs. |jSJ), ifTT]). and J22j) that 



-Z 



_ : 9Z 



dp \ i sps 
where (O)sps is defined by eq. (|2"o]) and 

1 9w fci 



k<i u k t 9/3 



We have used the relation 
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i<j 



(Bl) 



(B2) 



(B3) 



to go from the second line of eq. (|B1[) to the third. The 
derivative of Wki in the right hand side of eq. (|B2[) is 
calculated as 



dcu 



hi 



— h 

~W ~ dp I : 
= -<W,o (Vfci 



-^[Vfci+AVS] 



Sg klA [V kl +p 



dV k 



OB 



-pv kl 



kl 



Sg kh l[V k {+P 



dV, 



kl 



dp 



(B4) 



where E const and E are defined by eqs. (|2"T|) and (J35J), 
respectively. Equation ([24]) is obtained by substituting 
this equation into eq. (IB1[) . 



We next derive a formula for heat capacity (C) from 
cq. (|23|) . To this end, we first calculate the second deriva- 
tive of Z in the equation. We see from eq. (|B3[) that 



— TT 

9/3 2 -LJ- 



= E' 



— in 

9/3 i J-^ 



(B9) 



By taking the trace in the both hand slides of the equa- 
tion, we obtain 



9 2 Z /~, 2 dE> 
W~ \ ~9P 



(B10) 



SPS 



Equation (|25p is derived from this equation and eqs. 
(|Bip . and (|B8[) . Equation (f2T)| is readily obtained from 
eq. 
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